library(tidyverse)
library(janitor)
library(rsprite2)
library(scrutiny) 
library(tides)
library(purrr)
library(patchwork)
library(knitr)
library(kableExtra)

min_decimals <- function(x, digits = 2) {
  sprintf(paste0("%.", digits, "f"), x)
}

GRIM vs GRIM+GRIMMER vs GRIM+GRIMMER+Umbrella plot

Using {scrutiny} + {tides}

# dat_n_11_tides <-
#   expand_grid(mean = seq(from = 0.5, to = 7.5, by = 0.01),
#               sd = seq(from = 0, to = 4, by = 0.01)) |>
#   mutate(n_obs   = 11,
#          m_prec  = 2,
#          sd_prec = 2,
#          n_items = 1,
#          min_val = 1,
#          max_val = 7) |>
#   mutate(grim = pmap(list(as.character(mean), n_obs),
#                      grim)) |>
#   mutate(grimmer = pmap(list(as.character(mean), as.character(sd), n_obs),
#                         grimmer)) |>
#   mutate(tides = pmap(list(mean, sd, n_obs, min_val, max_val, n_items, sd_prec,
#                            calculate_min_sd = TRUE, verbose = FALSE),
#                       tides_single)) |>
#   unnest(grim) |>
#   unnest(grimmer) |>
#   unnest(tides)
# 
# dat_n_100_tides <-
#   expand_grid(mean = seq(from = 0.5, to = 7.5, by = 0.01),
#               sd = seq(from = 0, to = 4, by = 0.01)) |>
#   mutate(n_obs   = 100,
#          m_prec  = 2,
#          sd_prec = 2,
#          n_items = 1,
#          min_val = 1,
#          max_val = 7) |>
#   mutate(grim = pmap(list(as.character(mean), n_obs),
#                      grim)) |>
#   mutate(grimmer = pmap(list(as.character(mean), as.character(sd), n_obs),
#                         grimmer)) |>
#   mutate(tides = pmap(list(mean, sd, n_obs, min_val, max_val, n_items, sd_prec,
#                            calculate_min_sd = TRUE, verbose = FALSE),
#                       tides_single)) |>
#   unnest(grim) |>
#   unnest(grimmer) |>
#   unnest(tides)
# 
# write_rds(dat_n_11_tides, "results_grim_grimmer_tides_n_11_tides.rds")
# write_rds(dat_n_100_tides, "results_grim_grimmer_tides_n_100_tides.rds")

dat_n_11_tides <- read_rds("results_grim_grimmer_tides_n_11_tides.rds")
dat_n_100_tides <- read_rds("results_grim_grimmer_tides_n_100_tides.rds")

Plot and table functions

plot_helper <- function(dat, color_bounds = TRUE){
  
  dat <- dat |>
    mutate(label = case_when(mean < 1 | mean > 7 ~ "Impossible value flagged as possible",
                             sd < 0 | sd > 3.6 ~ "Impossible value flagged as possible", # Croucher (2004) loose heuristic of max SD as 60% of range (ie 6 * .6 = 3.6)
                             TRUE ~ "Possible value flagged as possible"))
  
  if(color_bounds == FALSE){
    
    # only GRIM-consistent values
    p_grim <- dat |>
      filter(grim) |>
      ggplot(aes(mean, sd)) +
      geom_point(shape = 15, size = 0.5, color = "grey35") + 
      theme_linedraw() +
      scale_y_continuous(breaks = scales::breaks_pretty(n = 8), 
                         limit = c(0, 4), 
                         expand = c(0.01, 0.01)) +
      scale_x_continuous(breaks = scales::breaks_pretty(n = 7), 
                         limit = c(0.5, 7.5), 
                         expand = c(0.01, 0.01)) +
      scale_color_viridis_d(begin = 0.3, end = 0.7) +
      ylab("Standard Deviation") +
      xlab("Mean") +
      ggtitle("GRIM consistent values") 
    
    # only GRIM and GRIMMER-consistent values
    p_grimmer <- dat |>
      filter(grim & grimmer) |>
      ggplot(aes(mean, sd)) +
      geom_point(shape = 15, size = 0.5, color = "grey35") +
      theme_linedraw() +
      scale_y_continuous(breaks = scales::breaks_pretty(n = 8), 
                         limit = c(0, 4), 
                         expand = c(0.01, 0.01)) +
      scale_x_continuous(breaks = scales::breaks_pretty(n = 7), 
                         limit = c(0.5, 7.5), 
                         expand = c(0.01, 0.01)) +
      scale_color_viridis_d(begin = 0.3, end = 0.7) +
      ylab("Standard Deviation") +
      xlab("Mean") +
      ggtitle("GRIM + GRIMMER consistent values") 
    
    # only GRIM and GRIMMER and TIDES-consistent values
    p_tides <- dat |>
      filter(grim & grimmer & tides_consistent) |>
      ggplot(aes(mean, sd)) +
      geom_point(shape = 15, size = 0.5, color = "grey35") +
      theme_linedraw() +
      scale_y_continuous(breaks = scales::breaks_pretty(n = 8), 
                         limit = c(0, 4), 
                         expand = c(0.01, 0.01)) +
      scale_x_continuous(breaks = scales::breaks_pretty(n = 7), 
                         limit = c(0.5, 7.5), 
                         expand = c(0.01, 0.01)) +
      ylab("Standard Deviation") +
      xlab("Mean") +
      ggtitle("GRIM + GRIMMER + TIDES consistent values")
    
    return(p_grim + 
             p_grimmer + 
             p_tides + 
             plot_layout(ncol = 1))
  }  
  
  if(color_bounds){
    
    # with coloring of impossible values
    # only GRIM-consistent values with colors
    p_temp_with_colors <- dat |>
      filter(grim) |>
      ggplot(aes(mean, sd, color = label)) +
      geom_point(shape = 15, size = 0.5) + # "grey20"
      theme_linedraw() +
      scale_y_continuous(breaks = scales::breaks_pretty(n = 8), 
                         limit = c(0, 4), 
                         expand = c(0.01, 0.01)) +
      scale_x_continuous(breaks = scales::breaks_pretty(n = 7), 
                         limit = c(0.5, 7.5), 
                         expand = c(0.01, 0.01)) +
      scale_color_viridis_d(begin = 0.3, end = 0.7) +
      ylab("Standard Deviation") +
      xlab("Mean") +
      ggtitle("GRIM consistent values") +
      theme(legend.position = "top",
            legend.direction = "vertical") +
      guides(color = guide_legend(override.aes = list(size = 4, ncol = 1), title = NULL))
    
    # Extract legend grob - magic chatGPT code no idea how it works
    legend_grob <- ggplotGrob(p_temp_with_colors)$grobs[[which(sapply(ggplotGrob(p_temp_with_colors)$grobs, function(x) x$name) == "guide-box")]]
    
    # Convert the legend grob into a patchwork-compatible object
    legend <- wrap_elements(grid::grobTree(legend_grob))
    
    # grim plot without legend
    p_grim_with_colors <- p_temp_with_colors +
      theme(legend.position = "none")
    
    # only GRIM and GRIMMER-consistent values with colors
    p_grimmer_with_colors <- dat |>
      filter(grim & grimmer) |>
      ggplot(aes(mean, sd, color = label)) +
      geom_point(shape = 15, size = 0.5) +
      theme_linedraw() +
      scale_y_continuous(breaks = scales::breaks_pretty(n = 8), 
                         limit = c(0, 4), 
                         expand = c(0.01, 0.01)) +
      scale_x_continuous(breaks = scales::breaks_pretty(n = 7), 
                         limit = c(0.5, 7.5), 
                         expand = c(0.01, 0.01)) +
      scale_color_viridis_d(begin = 0.3, end = 0.7) +
      ylab("Standard Deviation") +
      xlab("Mean") +
      ggtitle("GRIM + GRIMMER consistent values") +
      theme(legend.position = "none")
    
    # only GRIM and GRIMMER and TIDES-consistent values with colors
    p_tides_with_colors <- dat |>
      filter(grim & grimmer & tides_consistent) |>
      ggplot(aes(mean, sd, color = label)) +
      geom_point(shape = 15, size = 0.5) +
      theme_linedraw() +
      scale_y_continuous(breaks = scales::breaks_pretty(n = 8), 
                         limit = c(0, 4), 
                         expand = c(0.01, 0.01)) +
      scale_x_continuous(breaks = scales::breaks_pretty(n = 7), 
                         limit = c(0.5, 7.5), 
                         expand = c(0.01, 0.01)) +
      scale_color_viridis_d(begin = 0.3, end = 0.7, direction = -1) +
      ylab("Standard Deviation") +
      xlab("Mean") +
      ggtitle("GRIM + GRIMMER + TIDES consistent values") +
      theme(legend.position = "none")
    
    return(legend +
             p_grim_with_colors + 
             p_grimmer_with_colors + 
             p_tides_with_colors + 
             plot_layout(ncol = 1, heights = c(.1, 1, 1, 1)))
  }
}

# proportion of possible values, considering only those where mean is within the scale bounds and SD is within Croucher's (2004) loose heuristic of max SD as 60% of range (ie 6 * .6 = 3.6)
table_helper <- function(dat, min, max, max_sd){
  dat_feasible <- dat |>
    filter(mean >= min & mean <= max & sd >= 0 & sd <= max_sd)
  
  data.frame(check = c("GRIM", "GRIM+GRIMMER", "GRIM+GRIMMER+TIDES"),
             n_total = dat_feasible |>
               count(name = "n_total"),
             n_possible = c(dat_feasible |>
                              filter(grim) |>
                              count() |>
                              pull(n),
                            dat_feasible |>
                              filter(grim, grimmer) |>
                              count() |>
                              pull(n),
                            dat_feasible |>
                              filter(grim, grimmer, tides_consistent) |>
                              count() |>
                              pull(n))) |>
    mutate(proportion = janitor::round_half_up(n_possible/n_total, 3))
}

N=11

plot_helper(dat_n_11_tides, color_bounds = FALSE)

plot_helper(dat_n_11_tides, color_bounds = TRUE)

table_helper(dat_n_11_tides, min = 1, max = 7, max_sd = 6*.60) |>
  kable(align = "r") |>
  kable_classic(full_width = FALSE)
check n_total n_possible proportion
GRIM 46117 18085 0.392
GRIM+GRIMMER 46117 4139 0.090
GRIM+GRIMMER+TIDES 46117 3935 0.085

N=100

plot_helper(dat_n_100_tides, color_bounds = FALSE)

plot_helper(dat_n_100_tides, color_bounds = TRUE)

table_helper(dat_n_100_tides, min = 1, max = 7, max_sd = 6*.60) |>
  kable(align = "r") |>
  kable_classic(full_width = FALSE)
check n_total n_possible proportion
GRIM 141673 141673 1.000
GRIM+GRIMMER 141673 121201 0.855
GRIM+GRIMMER+TIDES 141673 112273 0.792

todo

  • wheres the code for the POMP plots? - in the /dev/relative folder